#---------------------------------------------------------
# Are Intermediary Constraints Priced?
# Du, Hébert, Huber
# The Review of Financial Studies 2022
#---------------------------------------------------------

#---------------------------------------------------------
# Bayesian model comparison as in Chib, Zeng, Zhao JF 2020 (CZZ)
# Based on functions in CZZ, as released in R package czzg
# Generate Tables 6, 7, A9
#---------------------------------------------------------

#---------------------------------------------------------
# Set-up
#---------------------------------------------------------

if (exists('script_path')) {
  rm(list = setdiff(ls(), 'script_path'))
} else {
  args = commandArgs(trailingOnly = T)
  if (length(args) > 0) {
    script_path <- args[1]  
    rm(list = setdiff(ls(), 'script_path'))
  } else {
    rm(list = ls())
    script_path <- '~/Dropbox/ForwardArbitrage/CodeReplication/'
  } 
}

# Inputs
source(paste0(script_path, 'RCode/ArbFunc.R'))
library(czzg) #download from <http://apps.olin.wustl.edu/faculty/chib/rpackages/>

load(paste0(script_path, 'OutputInterim/All_Rates_Dates_Basis_Returns.RData')) 
load(paste0(script_path, 'OutputInterim/Basis_and_Returns.RData')) 
load(paste0(script_path, 'OutputInterim/Factors_And_Test_Assets_Monthly.RData')) 

custom_czz <- function(xnames = xnames, wnames = NULL, data = data, trainpct = 0.1) 
{
  # Calculates log(marginal likelihood) as in CZZ 2020
  
  n = dim(data)[1]
  multiplier = (1 - trainpct)/trainpct
  nt = round(trainpct * n)
  n = n - nt
  kx = length(xnames)
  if (is.null(wnames)) {
    datnames = names(data) #to use this, must feed in data.FRAME
    wnames = setdiff(datnames, xnames)
  }
  kw = length(wnames)
  K = kx + kw
  X0 = as.matrix(data[1:nt, xnames, drop = FALSE])
  lampriorls = czzg::makelamprior(X0, K, kx, nt, multiplier)
  lamx0_ = lampriorls$lamx0_
  kappa = lampriorls$kappa
  data = data[-(1:nt), ]
  X = as.matrix(data[, xnames, drop = FALSE])
  if (kw > 0) {
    W = as.matrix(data[, wnames, drop = FALSE])
    out1 = gxmarg(X = X, K = K, lamx0_ = lamx0_, kappa = kappa)
    out2 = gwmarg(W = W, X = X, K = K)
    logmarg = out1 + out2
  } else {
    logmarg = gxmarg(X = X, K = K, lamx0_ = lamx0_, kappa = kappa)
  }
  return(logmarg)
}

#---------------------------------------------------------
# Data construction: monthly and quarterly series
#---------------------------------------------------------

## MONTHLY: HKM, market, int equity, basis, risk-free rate
## Note that the original risk-free rates are NOT lagged and are annualized, returns are from t-1 to t when taken from portfolio_monthly
## Carry returns are from t to t+1 if taken from basis_and_returns

hkm_monthly <- fread(paste0(script_path, 'DataSkeleton/AssetClasses/He_Kelly_Manela_Factors_monthly.csv')) #all rates in hkm are in UNITS
monthly_returns <- portfolio_monthly[, c('Date', 'USD_ois_1M', 'Tbill_1M', 'Mkt'), with = F] #all rates in portfolio_monthly are in UNITS
names(monthly_returns)[which(names(monthly_returns) == 'Mkt')] <- 'mkt_gross'

monthly_returns[, hkm_gross := hkm_monthly[, intermediary_value_weighted_investment_return][match(unlist(monthly_returns[, Date]), unlist(hkm_monthly[, yyyymm]))]]
monthly_returns[, hkm_factor := hkm_monthly[, intermediary_capital_risk_factor][match(unlist(monthly_returns[, Date]), unlist(hkm_monthly[, yyyymm]))]]
monthly_returns[, risk_free := shift(ifelse(is.na(USD_ois_1M), Tbill_1M, USD_ois_1M) / 12, type = 'lag')] 
monthly_returns <- monthly_returns[Date <= as.numeric(paste0(year(extended_end_date), month(extended_end_date)))]
monthly_returns[, `:=` (hkm_net = hkm_gross - risk_free, 
                        mkt_net = mkt_gross - risk_free)]
monthly_returns[, pc_net := portfolio_monthly[, pc1top6_ois_1M_fwd_3M][match(unlist(monthly_returns[, Date]), unlist(portfolio_monthly[, Date]))]]
monthly_returns[, carry_net := portfolio_monthly[, port_ois_ret_classic_1M_fwd_3M][match(unlist(monthly_returns[, Date]), unlist(portfolio_monthly[, Date]))]]

monthly_panel <- monthly_returns[complete.cases(monthly_returns[, c('carry_net', 'pc_net', 'hkm_net', 'mkt_net', 'hkm_factor')])]
monthly_panel <- monthly_panel[Date >= 201007, lapply(.SD, function(x) 100 * x), 
                               .SDcols = names(monthly_panel)[-1]] #state returns in percentage points to increase constant term

#---------------------------------------------------------
## QUARTERLY series: HKM / AEM are 3M returns ending in quarter-ends, basis returns / mkt are linear sum of three month-end returns (also t-1 to t since taken from portfolio_monthly)
## risk-free rates taken from portfolio_monthly are NOT lagged and are ANNUALIZED, but note that basis returns are already net

hkm_quarterly <- fread(paste0(script_path, 'DataSkeleton/AssetClasses/He_Kelly_Manela_Factors_quarterly.csv')) #all rates in hkm are in UNITS
names(hkm_quarterly) <- c('yyyyq', 'capital_ratio', 'hkm_factor', 'hkm_gross', 'lev_squared')
aem_quarterly <- fread(paste0(script_path, 'DataSkeleton/AssetClasses/aem_update.csv')) 
aem_quarterly[, yyyyq := as.numeric(paste0(substr(quarter, 1, 4), substr(quarter, 6, 99)))]
names(aem_quarterly)[which(names(aem_quarterly) == 'dlev')] <- 'aem_factor'

temp_quarterly <- portfolio_monthly[, c('Date', 'USD_ois_1M', 'Tbill_1M', 'Mkt', 'port_ois_ret_classic_1M_fwd_3M', 'pc1top6_ois_1M_fwd_3M'), with = F] #all rates in portfolio_monthly are in UNITS
temp_quarterly[, yyyyq := as.numeric(paste0(substr(Date, 1, 4), ifelse(substr(Date, 5, 6) %in% c('01', '02', '03'), 1,
                                                                       ifelse(substr(Date, 5, 6) %in% c('04', '05', '06'), 2,
                                                                              ifelse(substr(Date, 5, 6) %in% c('07', '08', '09'), 3,
                                                                                     ifelse(substr(Date, 5, 6) %in% c('10', '11', '12'), 4, 0))))))]
temp_quarterly[, `:=` (mkt_gross = sum(Mkt),
                       carry_net = sum(port_ois_ret_classic_1M_fwd_3M),
                       pc_net = sum(pc1top6_ois_1M_fwd_3M),
                       count_carry = sum(!is.na(port_ois_ret_classic_1M_fwd_3M)),
                       max_date = max(Date)), by = .(yyyyq)]
temp_quarterly <- temp_quarterly[count_carry == 3]
temp_quarterly <- temp_quarterly[Date == max_date]

quarterly_returns <- Reduce(function(x, y) merge(x, y, by = 'yyyyq'), list(temp_quarterly[, .(yyyyq, USD_ois_1M, Tbill_1M, mkt_gross, carry_net, pc_net)],
                                                                                  hkm_quarterly[, .(yyyyq, hkm_factor, hkm_gross)],
                                                                                  aem_quarterly[, .(yyyyq, aem_factor)]))
quarterly_returns[, risk_free := shift(ifelse(is.na(USD_ois_1M), Tbill_1M, USD_ois_1M) / 4, type = 'lag')] 
quarterly_returns <- quarterly_returns[yyyyq >= 20103 & yyyyq <= 20204]
quarterly_returns[, `:=` (hkm_net = hkm_gross - risk_free, 
                        mkt_net = mkt_gross - risk_free)]
quarterly_panel <- quarterly_returns[, lapply(.SD, function(x) 100 * x), .SDcols = names(quarterly_returns)[-1]] #state returns in percentage points to increase constant term

#---------------------------------------------------------
# Model comparison - Bayesian - with PC of near arbitrages
#---------------------------------------------------------

## (1) Monthly: market + int equity + cip vs. permutations
factor_names <- c('mkt_net', 'hkm_net', 'pc_net')
dat = data.frame(monthly_panel[, factor_names, with = F])
lm_mkt_int_pc <- vector()
name_mkt_int_pc <- vector()
for (l in 1:length(factor_names)) {
  x_names <- combn(factor_names, l)
  for (i in 1:ncol(x_names)) {
    lm_mkt_int_pc <- append(lm_mkt_int_pc, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
    name_mkt_int_pc <- append(name_mkt_int_pc, paste(x_names[, i], collapse = '_'))
  }
}
lm_vector <- lm_mkt_int_pc
lm_vector <- lm_vector - min(lm_vector)
post_mkt_int_pc <- sapply(lm_vector, function(x) exp(x) / sum(exp(lm_vector)))
names(post_mkt_int_pc) <- name_mkt_int_pc

## (2) Monthly: hkm factor + market + cip vs. permutations
factor_names <- c('mkt_net', 'hkm_factor', 'pc_net')
dat = data.frame(monthly_panel[, factor_names, with = F])
lm_mkt_hkm_pc <- vector()
name_mkt_hkm_pc <- vector()
x_names <- matrix(c('hkm_factor'), ncol = 1)
for (i in 1:ncol(x_names)) {
  lm_mkt_hkm_pc <- append(lm_mkt_hkm_pc, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_hkm_pc <- append(name_mkt_hkm_pc, paste(x_names[, i], collapse = '_'))
}
x_names <- matrix(c('hkm_factor', 'mkt_net', 'hkm_factor', 'pc_net'), ncol = 2) #HKM Factor is not traded so cannot be a test asset and hence must be a risk factor
for (i in 1:ncol(x_names)) {
  lm_mkt_hkm_pc <- append(lm_mkt_hkm_pc, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_hkm_pc <- append(name_mkt_hkm_pc, paste(x_names[, i], collapse = '_'))
}
x_names <- matrix(c('hkm_factor', 'mkt_net', 'pc_net'), ncol = 1)
for (i in 1:ncol(x_names)) {
  lm_mkt_hkm_pc <- append(lm_mkt_hkm_pc, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_hkm_pc <- append(name_mkt_hkm_pc, paste(x_names[, i], collapse = '_'))
}
lm_vector <- lm_mkt_hkm_pc
lm_vector <- lm_vector - min(lm_vector)
post_mkt_hkm_pc <- sapply(lm_vector, function(x) exp(x) / sum(exp(lm_vector)))
names(post_mkt_hkm_pc) <- name_mkt_hkm_pc

## (3) Quarterly: aem factor + market + cip vs. permutations
factor_names <- c('mkt_net', 'aem_factor', 'pc_net')
dat = data.frame(quarterly_panel[, factor_names, with = F])
lm_mkt_aem_pc <- vector()
name_mkt_aem_pc <- vector()
x_names <- matrix(c('aem_factor'), ncol = 1)
for (i in 1:ncol(x_names)) {
  lm_mkt_aem_pc <- append(lm_mkt_aem_pc, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_aem_pc <- append(name_mkt_aem_pc, paste(x_names[, i], collapse = '_'))
}
x_names <- matrix(c('aem_factor', 'mkt_net', 'aem_factor', 'pc_net'), ncol = 2) #AEM Factor is not traded so cannot be a test asset and hence must be a risk factor
for (i in 1:ncol(x_names)) {
  lm_mkt_aem_pc <- append(lm_mkt_aem_pc, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_aem_pc <- append(name_mkt_aem_pc, paste(x_names[, i], collapse = '_'))
}
x_names <- matrix(c('aem_factor', 'mkt_net', 'pc_net'), ncol = 1)
for (i in 1:ncol(x_names)) {
  lm_mkt_aem_pc <- append(lm_mkt_aem_pc, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_aem_pc <- append(name_mkt_aem_pc, paste(x_names[, i], collapse = '_'))
}
lm_vector <- lm_mkt_aem_pc
lm_vector <- lm_vector - min(lm_vector)
post_mkt_aem_pc <- sapply(lm_vector, function(x) exp(x) / sum(exp(lm_vector)))
names(post_mkt_aem_pc) <- name_mkt_aem_pc

#---------------------------------------------------------
# Model comparison - Bayesian - with Classic Carry
#---------------------------------------------------------

## (1) Monthly: market + int equity + cip vs. permutations
factor_names <- c('mkt_net', 'hkm_net', 'carry_net')
dat = data.frame(monthly_panel[, factor_names, with = F])
lm_mkt_int_carry <- vector()
name_mkt_int_carry <- vector()
for (l in 1:length(factor_names)) {
  x_names <- combn(factor_names, l)
  for (i in 1:ncol(x_names)) {
    lm_mkt_int_carry <- append(lm_mkt_int_carry, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
    name_mkt_int_carry <- append(name_mkt_int_carry, paste(x_names[, i], collapse = '_'))
  }
}
lm_vector <- lm_mkt_int_carry
lm_vector <- lm_vector - min(lm_vector)
post_mkt_int_carry <- sapply(lm_vector, function(x) exp(x) / sum(exp(lm_vector)))
names(post_mkt_int_carry) <- name_mkt_int_carry

## (2) Monthly: hkm factor + market + cip vs. permutations
factor_names <- c('mkt_net', 'hkm_factor', 'carry_net')
dat = data.frame(monthly_panel[, factor_names, with = F])
lm_mkt_hkm_carry <- vector()
name_mkt_hkm_carry <- vector()
x_names <- matrix(c('hkm_factor'), ncol = 1)
for (i in 1:ncol(x_names)) {
  lm_mkt_hkm_carry <- append(lm_mkt_hkm_carry, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_hkm_carry <- append(name_mkt_hkm_carry, paste(x_names[, i], collapse = '_'))
}
x_names <- matrix(c('hkm_factor', 'mkt_net', 'hkm_factor', 'carry_net'), ncol = 2) #HKM Factor is not traded so cannot be a test asset and hence must be a risk factor
for (i in 1:ncol(x_names)) {
  lm_mkt_hkm_carry <- append(lm_mkt_hkm_carry, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_hkm_carry <- append(name_mkt_hkm_carry, paste(x_names[, i], collapse = '_'))
}
x_names <- matrix(c('hkm_factor', 'mkt_net', 'carry_net'), ncol = 1)
for (i in 1:ncol(x_names)) {
  lm_mkt_hkm_carry <- append(lm_mkt_hkm_carry, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_hkm_carry <- append(name_mkt_hkm_carry, paste(x_names[, i], collapse = '_'))
}
lm_vector <- lm_mkt_hkm_carry
lm_vector <- lm_vector - min(lm_vector)
post_mkt_hkm_carry <- sapply(lm_vector, function(x) exp(x) / sum(exp(lm_vector)))
names(post_mkt_hkm_carry) <- name_mkt_hkm_carry

## (3) Quarterly: aem factor + market + cip vs. permutations
factor_names <- c('mkt_net', 'aem_factor', 'carry_net')
dat = data.frame(quarterly_panel[, factor_names, with = F])
lm_mkt_aem_carry <- vector()
name_mkt_aem_carry <- vector()
x_names <- matrix(c('aem_factor'), ncol = 1)
for (i in 1:ncol(x_names)) {
  lm_mkt_aem_carry <- append(lm_mkt_aem_carry, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_aem_carry <- append(name_mkt_aem_carry, paste(x_names[, i], collapse = '_'))
}
x_names <- matrix(c('aem_factor', 'mkt_net', 'aem_factor', 'carry_net'), ncol = 2) #AEM Factor is not traded so cannot be a test asset and hence must be a risk factor
for (i in 1:ncol(x_names)) {
  lm_mkt_aem_carry <- append(lm_mkt_aem_carry, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_aem_carry <- append(name_mkt_aem_carry, paste(x_names[, i], collapse = '_'))
}
x_names <- matrix(c('aem_factor', 'mkt_net', 'carry_net'), ncol = 1)
for (i in 1:ncol(x_names)) {
  lm_mkt_aem_carry <- append(lm_mkt_aem_carry, custom_czz(xnames = x_names[, i], wnames = NULL, data = dat, trainpct = 0.1))
  name_mkt_aem_carry <- append(name_mkt_aem_carry, paste(x_names[, i], collapse = '_'))
}
lm_vector <- lm_mkt_aem_carry
lm_vector <- lm_vector - min(lm_vector)
post_mkt_aem_carry <- sapply(lm_vector, function(x) exp(x) / sum(exp(lm_vector)))
names(post_mkt_aem_carry) <- name_mkt_aem_carry

#-----------------------------------------------------------
# Alpha test (GRS)
#-----------------------------------------------------------

monthly_pc_mean <- lm(pc_net ~ 1, data = monthly_panel)
quarterly_pc_mean <- lm(pc_net ~ 1, data = quarterly_panel)

monthly_mkt_return_unrest <- lm(pc_net ~ mkt_net, data = monthly_panel)
monthly_hkm_return_unrest <- lm(pc_net ~ hkm_net, data = monthly_panel)
monthly_mkt_hkm_return_unrest <- lm(pc_net ~ hkm_net + mkt_net, data = monthly_panel)
monthly_mkt_hkm_factor_unrest <- lm(pc_net ~ hkm_factor + mkt_net, data = monthly_panel)
quarterly_mkt_aem_factor_unrest <- lm(pc_net ~ aem_factor + mkt_net, data = quarterly_panel)

list_lm <- list(monthly_pc_mean, monthly_mkt_return_unrest, monthly_hkm_return_unrest, monthly_mkt_hkm_return_unrest, monthly_mkt_hkm_factor_unrest,
                quarterly_pc_mean, quarterly_mkt_aem_factor_unrest)
list_se <- list()
for (l in 1:5) {
  list_se[[l]] <- sqrt(diag(NeweyWest(list_lm[[l]], prewhite = F, lag = 12)))
}
for (l in 6:7) {
  list_se[[l]] <- sqrt(diag(NeweyWest(list_lm[[l]], prewhite = F, lag = 4)))
}

#-----------------------------------------------------------
# Output tables
#-----------------------------------------------------------
log_bayes <- file(paste0(script_path, 'OutputFinal/log_bayesian.txt'), open = 'wt')
sink(log_bayes, type = 'output', split = TRUE) 

## Print table preambles + panel A
library(stargazer)
stargazer(monthly_pc_mean, monthly_mkt_return_unrest, monthly_hkm_return_unrest, monthly_mkt_hkm_return_unrest, monthly_mkt_hkm_factor_unrest, 
          quarterly_pc_mean, quarterly_mkt_aem_factor_unrest,
          type = 'latex',
          se = list_se,
          title = 'Pricing Fwd CIP Returns with Intermediary Wealth',
          dep.var.caption = c('1M-fwd 3M classic pc forward CIP returns'),
          order = c(3, 4, 1, 2),
          covariate.labels = c('Market', 'Int. Equity', 'HKM Factor', 'AEM Factor'),
          column.labels = c('Monthly Returns', 'Quarterly Returns'),
          column.separate = c(5, 2),
          star.cutoffs = c(0.1, 0.05, 0.01),
          df = FALSE,
          label = '',
          omit.stat = c('adj.rsq', 'rsq', 'ser', 'f'),
          table.layout = "=lc#-t-s-",
          notes.label = 'Standard errors are in parentheses',
          notes.append = F,
          no.space = T,
          float = T)

## PANEL B - Bayesian table with PC
prob_dt <- data.table(`Factor Space` = c('{Market, Int. Equity,', 'Fwd. Ret. PC}', rep(NA, 5), '{HKM Factor, Market,', 'Fwd. Ret. PC}', rep(NA, 2),
                                         '{AEM Factor, Market,', 'Fwd. Ret. PC}', rep(NA, 2)),
                      Model = c('Market', 'Int. Equity', 'Market + Int. Equity', 'Fwd. Ret. PC', 'Market + Fwd. Ret. PC', 'Int.Equity + Fwd. Ret. PC',
                                'Market + Int. Equity + Fwd. Ret. PC', 
                                'HKM Factor', 'HKM Factor + Market', 'HKM Factor + Fwd. Ret. PC', 'HKM Factor + Market + Fwd. Ret. PC',
                                'AEM Factor', 'AEM Factor + Market', 'AEM Factor + Fwd. Ret. PC', 'AEM Factor + Market + Fwd. Ret. PC'),
                      Probability = c(post_mkt_int_pc['mkt_net'], post_mkt_int_pc['hkm_net'], post_mkt_int_pc['mkt_net_hkm_net'], post_mkt_int_pc['pc_net'], 
                                      post_mkt_int_pc['mkt_net_pc_net'], post_mkt_int_pc['hkm_net_pc_net'], post_mkt_int_pc['mkt_net_hkm_net_pc_net'], 
                                      post_mkt_hkm_pc['hkm_factor'], post_mkt_hkm_pc['hkm_factor_mkt_net'], post_mkt_hkm_pc['hkm_factor_pc_net'], post_mkt_hkm_pc['hkm_factor_mkt_net_pc_net'],
                                      post_mkt_aem_pc['aem_factor'], post_mkt_aem_pc['aem_factor_mkt_net'], post_mkt_aem_pc['aem_factor_pc_net'], post_mkt_aem_pc['aem_factor_mkt_net_pc_net']),
                      Subtotal = c(rep(NA, 2), sum(post_mkt_int_pc[(which(names(post_mkt_int_pc) %in% c('mkt_net', 'hkm_net', 'mkt_net_hkm_net')))]),
                                   rep(NA, 3), sum(post_mkt_int_pc[-(which(names(post_mkt_int_pc) %in% c('mkt_net', 'hkm_net', 'mkt_net_hkm_net')))]),
                                   rep(NA, 1), sum(post_mkt_hkm_pc[(which(names(post_mkt_hkm_pc) %in% c('hkm_factor_mkt_net', 'hkm_factor')))]),
                                   rep(NA, 1), sum(post_mkt_hkm_pc[-(which(names(post_mkt_hkm_pc) %in% c('hkm_factor_mkt_net', 'hkm_factor')))]),
                                   rep(NA, 1), sum(post_mkt_aem_pc[(which(names(post_mkt_aem_pc) %in% c('aem_factor_mkt_net', 'aem_factor')))]),
                                   rep(NA, 1), sum(post_mkt_aem_pc[-(which(names(post_mkt_aem_pc) %in% c('aem_factor_mkt_net', 'aem_factor')))])))
xtable_matrix <- xtable(prob_dt, caption = 'Bayesian Posterior Prob of SDF Models with Fwd. Ret. PC', label = 'bayesian_pc')
align(xtable_matrix) <- c('l', 'L{0.28\\textwidth}', 'L{0.45\\textwidth}', 'R{0.12\\textwidth}', 'R{0.11\\textwidth}') #c('l', rep('L{0.50\\textwidth}', 2), 'R{0.10\\textwidth}')
digits(xtable_matrix) <- 3
addtorow <- list()
addtorow$pos <- list(-1, 0, 3, 9, 13, nrow(prob_dt))
addtorow$command <- c(paste0('\\hline'),
                      paste0('\\hline'),
                      paste0('\\cdashline{2-4}'),
                      paste0('\\cdashline{2-4}'),
                      paste0('\\cdashline{2-4}'),
                      paste0('\\hline'))

print(xtable_matrix, include.rownames = FALSE, hline.after = c(-1, 7, 11, 15), add.to.row = addtorow,
      include.colnames = TRUE, floating = T, caption.placement = 'top') 

## Robustness - Bayesian table with Classic Carry
prob_dt <- data.table(`Factor Space` = c('{Market, Int. Equity,', 'Fwd. CIP Ret.}', rep(NA, 5), '{HKM Factor, Market,', 'Fwd. CIP Ret.}', rep(NA, 2),
                                         '{AEM Factor, Market,', 'Fwd. CIP Ret.}', rep(NA, 2)),
                      Model = c('Market', 'Int. Equity', 'Market + Int. Equity', 'Fwd. CIP Ret.', 'Market + Fwd. CIP Ret.', 'Int.Equity + Fwd. CIP Ret.',
                                'Market + Int. Equity + Fwd. CIP Ret.', 
                                'HKM Factor', 'HKM Factor + Market', 'HKM Factor + Fwd. CIP Ret.', 'HKM Factor + Market + Fwd. CIP Ret.',
                                'AEM Factor', 'AEM Factor + Market', 'AEM Factor + Fwd. CIP Ret.', 'AEM Factor + Market + Fwd. CIP Ret.'),
                      Probability = c(post_mkt_int_carry['mkt_net'], post_mkt_int_carry['hkm_net'], post_mkt_int_carry['mkt_net_hkm_net'], post_mkt_int_carry['carry_net'], 
                                      post_mkt_int_carry['mkt_net_carry_net'], post_mkt_int_carry['hkm_net_carry_net'], post_mkt_int_carry['mkt_net_hkm_net_carry_net'], 
                                      post_mkt_hkm_carry['hkm_factor'], post_mkt_hkm_carry['hkm_factor_mkt_net'], post_mkt_hkm_carry['hkm_factor_carry_net'], post_mkt_hkm_carry['hkm_factor_mkt_net_carry_net'],
                                      post_mkt_aem_carry['aem_factor'], post_mkt_aem_carry['aem_factor_mkt_net'], post_mkt_aem_carry['aem_factor_carry_net'], post_mkt_aem_carry['aem_factor_mkt_net_carry_net']),
                      Subtotal = c(rep(NA, 2), sum(post_mkt_int_carry[(which(names(post_mkt_int_carry) %in% c('mkt_net', 'hkm_net', 'mkt_net_hkm_net')))]),
                                   rep(NA, 3), sum(post_mkt_int_carry[-(which(names(post_mkt_int_carry) %in% c('mkt_net', 'hkm_net', 'mkt_net_hkm_net')))]),
                                   rep(NA, 1), sum(post_mkt_hkm_carry[(which(names(post_mkt_hkm_carry) %in% c('hkm_factor_mkt_net', 'hkm_factor')))]),
                                   rep(NA, 1), sum(post_mkt_hkm_carry[-(which(names(post_mkt_hkm_carry) %in% c('hkm_factor_mkt_net', 'hkm_factor')))]),
                                   rep(NA, 1), sum(post_mkt_aem_carry[(which(names(post_mkt_aem_carry) %in% c('aem_factor_mkt_net', 'aem_factor')))]),
                                   rep(NA, 1), sum(post_mkt_aem_carry[-(which(names(post_mkt_aem_carry) %in% c('aem_factor_mkt_net', 'aem_factor')))])))
xtable_matrix <- xtable(prob_dt, caption = 'Bayesian Posterior Prob of SDF Models with Fwd. CIP Return', label = 'bayesian_carry')
align(xtable_matrix) <- c('l', 'L{0.28\\textwidth}', 'L{0.45\\textwidth}', 'R{0.12\\textwidth}', 'R{0.11\\textwidth}') #c('l', rep('L{0.50\\textwidth}', 2), 'R{0.10\\textwidth}')
digits(xtable_matrix) <- 3
addtorow <- list()
addtorow$pos <- list(-1, 0, 3, 9, 13, nrow(prob_dt))
addtorow$command <- c(paste0('\\hline'),
                      paste0('\\hline'),
                      paste0('\\cdashline{2-4}'),
                      paste0('\\cdashline{2-4}'),
                      paste0('\\cdashline{2-4}'),
                      paste0('\\hline'))

print(xtable_matrix, include.rownames = FALSE, hline.after = c(-1, 7, 11, 15), add.to.row = addtorow,
      include.colnames = TRUE, floating = T, caption.placement = 'top')

closeAllConnections()
file.show(paste0(script_path, 'OutputFinal/log_bayesian.txt'))

